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ABSTRACT 

We model the generation of a magnetic field in a protostellar disc using an a- 
dynamo and perform axisymmetric magnetohydrodynamics (MHD) simulations of a 
T Tauri star. We find that for small values of the dimensionless dynamo parameter 
otd the poloidal field grows exponentially at a rate a oc flxy/cui, before saturating 
to a value oc ,/cq). The dynamo excites dipole and octupole modes, but quadrupole 
modes are suppressed, because of the symmetries of the seed field. Inital seed fields 
too weak to launch MHD outflows are found to grow sufficiently to launch winds with 
observationally relevant mass fluxes of order 10 _9 MQ/yr for T Tauri stars. For large 
values of ad magnetic loops are generated over the entire disc. These quickly come 
to dominate the disc dynamics and cause the disc to break up due to the magnetic 
pressure. 

Key words: accretion, accretion discs - dynamo - MHD - stars: magnetic fields - 
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1 INTRODUCTION 

Magnetic fields play a key role in transporting angular mo¬ 
mentum away from protostellar systems, thereby allowing 
matter to accrete on time scales consistent with models of 
protostellar evolution. Turbulent viscosity in the accretion 
disc, most likely generated by the magnetorotational insta¬ 
bility (MRI), plays a key role in setting the accretion rate 
in the disc (Balbus & Hawley 1991). Astrophysical jets and 
winds are the dominant mechanisms for transporting angu¬ 
lar momentum away from protostellar systems. Magnetic 
fields play a central role in the launching and collimation 
of these outflows. The origin of these large-scale, ordered, 
magnetic fields required for sustaining jets and winds is at 
present not well understood. 

One possibility is that the magnetic field is advected 
inwards from the interstellar medium or ambipolar diffusion 
(see Shu, Adams & Lizano 1987 for review). Another pos¬ 
sibility is that the field is generated locally via a dynamo 
mechanism. It is interesting therefore, to simulate the ef¬ 
fects of a magnetic disc dynamo on protostellar evolution. 

Different theoretical dynamo models have been pro¬ 
posed. The mean field dynamo (Moffatt, 1978) relies on a 
magnetic field that is disordered on small scales interacting 
through the mean motion of the disc on large scales. Tout & 


Pringle’s (1992) mechanism relied on the Parker instability, 
MRI and magnetic reconnections. 

Observations suggest that a dynamo mechanism is re¬ 
sponsible, at least in part, for the magnetic fields near pro¬ 
toplanetary stars. Empirical evidence suggests there is a 
correlation between the magnitude of large scale magnetic 
fields and seed fields, suggesting a related dynamo mech¬ 
anism is responsible for generating them (Vidotto et al, 
2014). Rotational modulation of Zeeman signature in GQ 
Lup (Donati et al. 2012) and DN Tau (Donati et al. 2013) 
suggests that the magnetic fields of T Tauri stars are gen¬ 
erated by non-stationary dynamos. 

Much work has been done to understand the dynamo 
process in the context of generating turbulent magnetic 
fields to model the MRI. Stratified shearing box simulations 
with zero net vertical magnetic flux exhibit a flipping of the 
mean toroidal magnetic field flux on time scales of roughly 
10 orbits (Brandenburg et al. 1995; Shi et al. 2010; Davis 
et al. 2010, Simon et al. 2012). This effect is found to be 
more sporadic as the amount of net vertical magnetic flux 
is increased (Bai & Stone 2013). Full 3D simulations (Stone 
et al. 1996) of turbulent field evolution largely agreed with 
these shearing box models. Del Zanna et al. (2014) imple¬ 
mented a disc dynamo in the GRMHD code ECHO and 
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found exponential growth of toroidal magnetic fields from 
a seed poloidal field. Sadowski et al. (2014) used a sub-grid 
model in the GRMHD code KORAL to simulate a dynamo 
in a thick disc in 2D and found accretion rates were compa¬ 
rable to 3D simulations. The disc dynamo played a key role, 
by allowing for the regeneration of poloidal magnetic field 
which drove the disc to a state characteristic for a 3D MRI 
disc. Schober et al (2012) have shown that dynamo mecha¬ 
nisms can amplify magnetic fields during the formation of 
primordial halos by converting kinetic energy from turbu¬ 
lence into magnetic energy, von Rekowski (2003) found that 
an a 2 dynamo can launch outflows even without an exter¬ 
nal field. Bardou et al. (2001) simulated the back-reaction 
of outflows on the disc dynamo and found a vertical ve¬ 
locity can enhance the action of the dynamo. In a series 
of papers Torkelsson & Brandenburg (1994a; 1994b; 1994c) 
studied the parity of generated fields and growth rates of 
dipolar and quadrupolar magnetic fields from linear and 
non-linear cc-dynamo effects. They found their growth to 
be chaotic for large enough dynamo number. Stepanovs, 
Fednt & Sheikhnezami (2014) generated episodic outflows 
with a time dependent mean field disc dynamo. 

In this paper we perform axisymmetric MHD simula¬ 
tions of an a-dynamo operating inside a protostellar accre¬ 
tion disc. We provide an initial weak poloidal seed magnetic 
field and track its growth and the corresponding MHD out¬ 
flows. Our goal is to undertand how the disc dynamo affects 
the corresponding outflows and the late time magnetic field 
structure. The structure of this paper is as follows. In Sec¬ 
tion 2 we review the disc dynamo contribution to the equa¬ 
tions of motion and describe our MHD code. In section 3 
we present results, describing how outflows and magnetic 
field structure around the star are affected by the initial 
seed field, the dynamo parameter ad, the disc viscosity a v 
and the magnetic diffusivity a v - We conclude in Section 4 
and discuss implications for future work and observational 
prospects for classical T Tauri stars (CTTS). 


2 DYNAMO MODEL 
2.1 Basic Equations 

The plasma flows are assumed to be described by the equa¬ 
tions of non-relativistic magnetohydrodynamics (MHD). In 
a non-rotating reference frame the equations are 

g+V-(pv) = 0, (la) 

^ + V ■ T = ps, (lb) 

®-Vx(vxB) + Vx (»jtV x B) + V x (AB) = 0 , (lc) 

+ V ■ (pvS 1 ) = Q , (Id) 

where p is the mass density, S is the specific entropy, v is the 
flow velocity, B is the magnetic field, T is the momentum 
flux density tensor, E = — vxB/c+J)(VxB/c is the electric 
field, Q is the rate of change of entropy per unit volume due 


to heating, rjt is the turbulent magnetic diffusivity and c is 
the speed of light. 

The dynamo effect is generated by the turbulent fluc¬ 
tuations of the velocity 5v. Above 

A = -^(<5v ■ V x <5v) = Zn K a d , (2) 

is the product of the decorrelation time of the velocity fluc¬ 
tuations t and the mean lielicity of the turbulence. We 
parametrize this by an a prescription, with the strength 
of the effect parametrized by the a-dynamo parameter ad- 
See Kulsrud (1999) for review of the mean field theory dy¬ 
namo. ad is a dimensionless number that we vary in our 
simulations and the prefactor of the height in the disc Z 
times the local Keplerian angular velocity f Ik was chosen 
by dimensional analysis. 

We assume that the heating is offset by radiative cool¬ 
ing so that Q — 0. Also, g = — (GM/r 2 ) f is the gravita¬ 
tional acceleration due to the central mass M. We model 
the fluid as a non-relativistic ideal gas with equation of state 


5 = In 



( 3 ) 


where p is the pressure and 7 = 5/3. 

We use cylindrical coordinates ( R , </, Z ) as these are the 
coordinates used by our code. Sometimes it will be useful 
to also refer to spherical coordinates (r, 9, (j> ) as well. 

We consider that both the viscosity and the magnetic 
diffusivity of the disc plasma are due to turbulent fluctu¬ 
ations of the velocity and magnetic field. We assume that 
the molecular transport coefficients can be replaced by the 
turbulent coefficients. The turbulent coefficients are param¬ 
eterized using the a-model of Shakura and Sunyaev (1973). 
The turbulent kinematic viscosity is 

vt = a„^, (4) 

11k 

where c s is the midplane sound speed and q„ ^ 1 is a 
dimensionless constant. Similarly, the turbulent magnetic 
diffusivity is 

Vt = Oi v -y- , (5) 

\Lk 

where a v is another dimensionless constant. The ratio, 


■p = — 
Ot-r) 


( 6 ) 


is the magnetic Prandtl number of the turbulence in the 
disc. Shearing box simulations of MRI driven MHD turbu¬ 
lence in discs indicate that V ~ 1 (Guan & Gammie 2009). 
We take a v = a v = 0.1 unless otherwise stated. 

We define the dynamo number 

Nd = AnKtf^ad (7) 

fit al 


Dynamo effects are thought to be important for Nd > N c 
for some critical dynamo number (Stepinski & Levy 1988). 
The momentum flux density tensor is given by 

'T x _i_ , /B 2 BiB k \ , . 

Tik = POik + pViVk + I —Oik - I + Tik , (8) 


where Tik is the viscous stress contribution from the turbu¬ 
lent fluctuations of the velocity. As mentioned, we assume 
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that these can be represented in the same way as the col- 
lisional viscosity by substitution of the turbulent viscosity. 
The leading order contributions to the viscous stress arise 
from large velocity gradients. In a Keplerian type disc these 
are dominated by the azimuthal terms v$ ~ Vk- However, 
we include all viscous terms involving vr and v^. The lead¬ 
ing order contribution to the momentum flux density from 
turbulence are therefore 


= —v t pR 


dQ 


tr$ = 

TZR = —Vtp 

T 4>$ = 


= -2 v t p 


dR ’ 
dv R 

Fz ’ 

Vr 


Tz<t> = 

Trr = 


—vtpR 
-2 v t p 


dQ 

dz 

dv R 


R 


(9) 


where = v^/R is the angular velocity of the gas. 

The transition from the viscous disc to the non-viscous 
corona is handled by multiplying the viscosity by a dimen¬ 
sionless factor £(p) which varies smoothly from £ = 1 for 
p ^ Pd = 0.75p(.R, Z = 0) to ^ = 0 for p ^ 0.25p,j as de¬ 
scribed in Appendix B of Lii, Romanova, & Lovelace (2012). 


2.2 Disc Averaged Equations 


To gain some insight into the dynamo process, consider a 
thin, axisymmetric, Keplerian disc with diffusivity p and 
radial infall velocity vr operating an a-dynamo. We will 
ignore any vertical motion vz and neglect the Z dependence 
of ff and vr. Breaking the induction equation (JTc]) into (R, 
cf), Z) components, we find 


1 a 

(»\ 

d 

RdZ 

[dt) 

~dZ 


— AB ^ + 


Vr <93> 

ADlm 


1 (nJL (ld*\ d 2 3>\ 

R\ dR\RdR) dZ 2 ) ’ 

(10a) 


dt. 


tmd^_ d_(_Ad^ dBfi\ 
dRdZ + dZ\ RdZ +V dZ ) 


d_ 

dR 


( A&$_ 
\RdR 


+ vrB$ 


V d{RBD\ 
R dR ) ’ 


(10b) 


9 

(d*\ 

d 

dR 

[dt] 

~dR 


d 3< 

AB^R - v R-qr 


+ V 



\R~dR) + ~dZ 2 ) ’ 


(10c) 


where 'I' is the magnetic flux function. Two types of solu¬ 
tions are possible, depending on the symmetries of 3' and 
B,/, - a “dipole” like solution where 3>(Z) = 'i(-Z) and 
B,f 1 {Z) = —B^—Z) and a “quadrupole” solution where 
3>(Z) = -'F(-Z) and B${Z) = B^-Z) (see Kulsrud 
2005). To solve the problem exactly we make a gauge choice, 
for example 31(0) = 0 and boundary conditions at Z = izh 
such as 93t fdZ = 0. Intuitively then, one may think of the 
flux function as a cosine in the dipole like case and as a sine 
in the quadrupole case 

* ( z)-{* 0COB ,(f) “ dipole ”’ (11) 

( 3f 0 sin (^) “quadrupole”. 


We work in the thin disc, h <C R, approximation. The 
variation in the flux function over the half disc thickness 
is therefore in each case A 31 ss 31o but the variation in 
its derivative for each case is A3//AZdi p ss 2A3//A^ qua d. 
Therefore integrating over the upper half part of the disc, 
and replacing all quantities with their vertically averaged 
quantities J Q h B^dz = hBfi 1 —^§§= IiBr and ignoring 
terms involving Bz, the R and <j> components give us the 
coupled system 

^ = (|Cr + C 2 a d )nB R - A_B^ (12a) 


= C 2 a d ^B^ - AL Br . (12b) 

where C\ = 1(—1) and C 2 = 2(—1) in th e d ipole 
(quadrupole) case. Using r\ = a v flh 2 , we combine (121 and 
find 

* S «+2 


dt 2 ' dt 

+ fir] — C 2 ot d (3/2Ci + C 2 otd) +] 12 Br = 0. 
Taking Br oc exp(<rf2t) we find growth rates 

a = —otn ± \JC 2 a d (3/2Ci + C 2 a d ). 
Taking a d -C 1 we approximate 

I ±\/3 a d “dipole”, 


a + a „ 


| ±^/3a d /2 “quadrupole” 


Substituting back into (121 we find 



A^Bj, “dipole”, 
^^B^ “quadrupole”. 


(13) 


(14) 


(15) 


(16) 


The toroidal field component is primarily generated by the 
differential rotation of the disc. It is also the dominant con¬ 
tribution to the magnetic field. Therefore, at late times 
when the matter and magnetic pressure in the disc are in 
equilibrium, we expect the toroidal component to saturate. 
This would imply that the radial field should also saturate, 
to a value proportional to yjotd- For this reason, we scale 
all our magnetic field plots to fictd . Also, for ad,ot v <S 1 
the growth rate a oc fiod and we scale all times by a fac- 
tor of a d ' . We note that requiring a positive growth rate 
a > 0 implies N d > 1/3. Though we have neglected order 
one numbers, this provides an inuitive picture for why there 
is a critical dynamo number - if the diffusion time scale is 
too fast compared to the dynamo time scale then the field 
will diffuse away and never grow. Stepinski & Levi (1988) 
found purely growing dynamo modes for dynamo number 
R m = ( R/h ) 3 ' 2 a d /a 2 ~ 100 which agrees with this result 
up to some 0( 1) numbers. 

Above we neglected terms that depended on d/dR to 
make the exponential growth of the field via the dynamo 
apparent. However, these terms are important for deter¬ 
mining the effects of viscosity and diffusivity in the disc. In 
particular, taking vr = v/R and v,p constant and expand¬ 
ing (101 will yield advection like terms oc (—v + ^d'it/dR 
and (—v — r/)dB^/dR. rj and v may cooperate or compete in 
transporting different field components. We show this em¬ 
pirically in our simulations where the effects of the dynamo 
can be changed by varying the Prandtl number V = a„/a,. 
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2.3 Initial Conditions 

2.3.1 Magnetic Field & Dynamo 


Our simulations use a disc field (Zanni, 2007) as a seed field 
for the dynamo, defined by the poloidal flux function 


I'disc = ^B 0 Rl 


RoJ 


3/4 


,5/4 


(m 2 + Z 2 /R 2 ) 5/8 ’ 


(17) 


where the poloidal magnetic field components can be com¬ 
puted via 


B z 


1 d'H 
RdR’ 


Br = - 


1 d'H 
RdZ' 


(18) 


The parameter m determines how much the field lines bend 
in the R-Z plane with the limit m —► oo corresponding to 
purely vertical field lines. In this study we set m = 1. 

The field generated by the dynamo can be expanded 
using a multipole expansion 


U/ — /4fl i p 4/di p T //quad 4/quad T //oct 4/ oc t “t“ (10) 


where 


^dip = 


R 2 

(R 2 + Z 2 )l’ 


is the dipole contribution to the potential, 

_ 3 R 2 Z 
^ ~ 4 (^2 + ^ 2 ) 5/2 ’ 

is the quadrupole contribution and 

_R 2 (4Z 2 -R 2 ) 
oct ~ 2 (R 2 + Z 2 y/ 2 ’ 


( 20 ) 


( 21 ) 


( 22 ) 


is the octupole contribution. These are expected to be the 
dominant modes excited by the disc dynamo. The disc dy¬ 
namo is given by 


A(R, Z) 


0 

a d zJ^f e - (pd y* /p)2 tanh( 


R -Rdyn \ 
AH ) 


R Rdyn , 

R ' Rdyn • 
(23) 


In our simulations we have set Rd yn = 5, pdyn = 0.1 and 
A R = 1. This is the form of the a—dynamo described in 
|2| but with an exponential factor to continuously cut it off 
from operating in the corona and a hyperbolic tangent to 
continuously cut it off at an inner radius Rdyn- We found it 
necessary to impose an inner disc cutoff to avoid numerical 
instabilities at the surface of the star. 


2.3.2 Matter Distribution 


Initially the matter of the disc and corona are assumed to 
be in mechanical equilibrium (Romanova et al. 2002). The 
initial density distribution is taken to be barotropic with 


p(p) 


p/T diac p> p b and R ^ R d , 
p/T COI p < p b or R 4 . R d , 


(24) 


where pb is the level surface of pressure that separates the 
cold matter of the disc from the hot matter of the corona 
and R d is the inital inner disc radius. At this surface the 
density has an initial step discontinuity from value p/Tdi S c 
to p/T cor . 



R 


Figure 1. Sparse version of grid used in the simulations, showing 
every 11th gridline in the R direction and every 10th gridline in 
the Z direction. 


Because the density distribution is barotropic, the ini¬ 
tial angular velocity is a constant on coaxial cylindrical sur¬ 
faces about the Z— axis. Consequently, the pressure can be 
determined from the Bernoulli equation 


F(p) + 4> + $ c = Bo , (25) 

here Bo is Bernoulli’s constant,4? = — GM/\/R 2 + Z 2 is the 
gravitational potential, with GM = 1 in the code, <E> C = 
rdr cu 2 (r) is the centrifugal potential, and 


F(p) 


Tdisc In (p/p b ) p> Pb and R^ R d , 

Tcor In (p/p b ) p < Pb or R^ R d . 


(26) 


We take the Bernoulli parameter Bo = 3 x 10 -4 in our 
simulations. The initial half thickness of the disc h/R = 0.1 
at the initial inner disc edge R d = 5. 


2.3.3 Angular Velocity 

The initial angular velocity of the disc is slightly sub- 
Keplerian, 

12 = (1 - 0.003)^(77) R>R d , (27) 

Inside of R d , the matter rotates rigidly with angular velocity 

12 = (1 — 0.003)12 k (Rd) R^R d . (28) 

The corotation radius R c is the radius where the angular 
velocity of the disc equals that of the star; that is, R c = 
(GM/Q 2 ) 1 / 3 . In this study we have chosen R c = 2. 


2.4 Boundary Conditions 

Our simulation region has three boundaries: the axis, the 
surface of the star and the external boundaries. For each 
dynamical variable we impose a boundary condition consis¬ 
tent with our physical assumptions. 
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Parameters 

Symbol 

Value 

mass 

Mo 

1.59 x 10 33 g 

length 

Ro 

1.50 x 10 12 cm 

magnetic field 

Bo 

8.04 x 10 -2 G 

time 

Po 

3.53 x lO" 2 y 

velocity 

VO 

8.42 x 10 7 cm/s 

density 

PO 

7.24 X 10~ 15 g/cm 3 

accretion rate 

M 0 

2.72 X 10“ 7 M 0 /yr 

temperature 

To 

4.26 X 10 7 K 

dipole strength 

M dip 

2.69 x 10 35 G cm 3 

quadrupole strength 

Mquad 

4.04 x 10 47 G cm 4 

octupole strength 

/^oct 

6.05 x 10 59 Gem 5 


Table 1. Mass, length, and magnetic field scales of interest and 
the corresponding scales of other derived quantities for a CTTS. 
One can obtain the physical values from the simulation values by 
multiplying by the corresponding dimensional quantity above. 

We assume axisymmetry about the axis. On the star 
and the external boundaries we want to allow fluxes and 
impose free boundary conditions dJ-/dn = 0 where J- is 
a dynamical variable and n is the vector normal to the 
boundary. 

At the external boundary along the edge of the disc 
we allow new matter to flow into the simulation region. We 
impose the condition that the matter must be accreting 
v r < 0. In the coronal region, we prescribe outflow condi¬ 
tions and allow matter, entropy and magnetic flux to exit 
the simulation region. 

The star is taken to be cylindrical in shape with radius 
R„ = 1 and height Z* = 2. Free boundary conditions are 
imposed on the variables p, p and B. In addition, we require 
that matter on the stellar boundary flow into the star. We 
treat the corner of the star by averaging over the nearest 
neighbour cells in the R and Z directions. 

Our simulation uses a grid Nr x Nz = 154 x 230 cells. 
The star has a radius of 1 in units of the simulation and is 
cylindrical in shape. It extends 10 grid cells above and be¬ 
low the equatorial plane. In the R-direction, the first 60 grid 
cells have length dR = 0.05. Later cell lengths are given re¬ 
cursively by dRi +1 = 1.025<iR; . Similarly, in the Z-direction 
the first 30 grid cells above and below the equatorial plane 
have length dZ = 0.05. Later cell lengths are given recur¬ 
sively by dZj +1 = 1.025 dZj. We show a sparse version of 
our grid, showing every 11th gridline in the R direction and 
every 10th gridline in the Z direction in Fig. [l] 

2.5 Dimensional Variables 

The MHD equations are written in dimensionless form so 
that the simulation results can be applied to different types 
of stars. We assume that the central object is a CTTS with 
mass M* = 0.8 Mg, a radius R* = 2 Rq and a magnetic 
field with magnitude B* = 3 x 10 3 G on the stellar surface, 
which is typical for the magnitude of the stellar dipole. We 
define dimensionful quantities with a 0 subscript, to denote 
typical values of physical parameters at a reference radius 
Ro. The mass scale is set by the stellar mass so Mo = M*. 
The reference length, Rq = 0.1 AU, is taken to be the scale 


of a typical inner disc radius. Assuming a stellar dipole field, 
the magnetic field strength Bo = -B,(R„/Ro) 3 - The mass, 
length and magnetic field scales allow us to define all other 
dimensionful quantities. 

The reference value for the velocity is the Keplerian ve¬ 
locity at the radius Ro, Vo = (GMo/Ro) 1 ^ 2 ■ The reference 
time-scale is the period of rotation at Ro, Po = 2-kRq/vo- 
From the relation pov p = Bo/47r, we define the reference 
density po of the disc. This also defines a reference pres¬ 
sure The reference mass accretion rate is Mo = 4-7rpoVoRo- 
The reference temperature To = Vq/2 x mH/ks where rriH 
is the atomic mass of hydrogen and fes, the Boltzmann 
constant, is obtained by taking the ratio of the reference 
pressure po = Bq / 8tt and the reference density. The dimen¬ 
sions of the coefficients in the multipole expansion of the 
flux function are obtained by scaling the reference field by 
the appropriate power of R 0 - p d i P = B 0 Rq, Mquad = PdipRo 
and Poet — Pdip Ro * 

Results obtained in dimensionless form can be applied 
to objects with widely different sizes and masses. However, 
the present work focuses on CTTS with the typical values 
shown in Table [l] One can obtain dimensionful quantities 
from simulation results by multiplying by the appropriate 
dimensionful quantity above. A more detailed description of 
our code and numerical methods can be found in Koldoba 
et al. (2015). 


3 RESULTS 

The a-dynamo converts a poloidal magnetic field into a 
toroidal magnetic field and vice versa. We provide an initial 
poloidal seed field. The differential disc rotation causes a 
twisting of this poloidal field, producing a toroidal field. The 
a dynamo in turn converts this toroidal field into poloidal 
field. This is a positive feedback loop and the poloidal 
field continues to grow until other physical effects halt the 
growth. This could have interesting observational implica¬ 
tions, as seed fields too weak to launch MHD outflows can 
grow to values which can support observable outflows. 

We find two qualitatively different regimes in our dy¬ 
namo simulations - a weak dynamo regime where dynamo 
effects are only important in the inner disc and a strong 
dynamo regime where dynamo effects affect the entire disc 
and the magnetic field quickly comes to dominate the disc 
dynamics. We discuss our results for each of these cases 
below. 


3.1 Weak Dynamo 

Below we show poloidal plots of disc density p, magnetic 
flux contours T and poloidal velocity vectors v p of a typ¬ 
ical run in the weak dynamo regime (Fig. [2|. We see that 
magnetic flux of the disc increases as a function of time, 
and the flux builds up on the star as it advects inwards 
with the accreting matter. The dynamo only changes the 
field geometry in the inner part of the disc R < 10 where 
loops form, even though the dynamo is operating for all 
radii R > Rd yn = 5. In the outer parts of the disc field lines 
become more inclined relative to their initial configuration, 
just like in the case where the dynamo is not operating. 
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Figure 2. Density p (color), poloidal field lines 'I' (white) and poloidal velocity r p for the case lV d = 10. t = 0 shows the initial large 
scale field. By t = 100 the dynamo has changed the field structure in the inner part of the disc. By t = 200 the field has developed a 
parallel dipole structure near the star, which interacts with the large scale field to produce an X-point at R rj 6 . By t = 300 the dipole 
structure has pushed out the inner part of the disc. 


3.1.1 Stellar Multipoles 

Any magnetic field can be expanded in a general multi¬ 
pole expansion. As our code is in cylindrical coordinates, 
projecting out the multipole components is slightly more 
involved than taking the inner product of the flux function 
with the Legendre polynomials. We show how to extract 
the dipole /i dip , quadrupole /r quad and octupole p oct mo¬ 
ments from our simulations in Appendices |A1| - |A3| We 
use this method to study the time dependence of the mul- 
tipole moments. This gives us a sense of the geometry of 
the generated magnetic field around the star. This is im¬ 
portant because measurements using the Zeeman-Doppler 
technique have shown that magnetic fields around young 
stars are often complex, exhibiting dipole and higher order 
fields (e.g Donati & Collier Cameron 1997; Donati et al. 
1999; Jardine et al. 2002). 


Runs were performed with seed field amplitudes rang¬ 
ing between Bo = 10~ 2 — 10°. These correspond to mid¬ 
plane plasma (3 = 10 1 — 10 5 . For f) 2> 10 1 , corresponding 
to seed fields Bo <C 1, the results were largely independent 
of the initial magnetic field. For larger field values, the dy¬ 
namo had comparatively little effect, and the simulation 
behaved as for the case a d = 0 where the disc dynamo is 
explicitly turned off. Growth of the dipole and higher order 
moments around the star is due to advection and rearrange¬ 
ment of the field, and not the generation of new magnetic 
flux. Cowling’s antidynamo theorem (1934) states that a 
self-sustained, axisymmetric field cannot be maintained - in 
essence non axisymmetric velocity components must coun¬ 
teract the diffusion of the field. Our simulations run on the 
order of 600 orbits at R = 1, corresponding to roughly 6 
orbits at R = 20. The typical diffusion time t ~ R 2 /r/t ^ 1, 
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meaning we can evade the antidynamo theorem because we 
have not waited sufficiently long for the field to diffuse away. 
In fact, after t = 600 the total magnetic flux threading the 
star and disc decreased roughly 5% from its initial value. 
The flux on the star increased roughly 100% and the flux 
threading the disc decreased approximately 20%. 

For fixed diffusivity a v = 0.1 and seed field amplitude 
B 0 = 0.1 we plot the dipole (Fig. |3|, quadrupole (Fig. [4j) 
and octupole moments (Fig. [5| as a function of time for 
dynamo numbers Nd = 1,5,10 (ad = 0.01,0.05 and 0.1). 
We have scaled the time by a d and the dipole moment 


by a 


- 1/2 


as suggested by the linear analysis of Section 


2.2 


Empirically, we see this scaling yields similar behaviour for 
runs with different dynamo strength parameter when ad -C 
1. We find qualitatively different behaviour of the dipole, 
quadrupole and octupole moments. The dipole moment 
grows to a value of approximately /idip ~ —iV d ' . It sat¬ 
urates to this value on time scales r ~ a^ 1 ' /2 fl.K'(Rd y n) _1 . 
The quadrupole moment oscillates about /i qua d = 0, which 
is what we would expect given the symmetry of the ini¬ 
tial seed field. The magnitude of these oscillations increases 
with increasing dynamo number. At late times the ampli¬ 
tude of the oscillations dies down. The octupole moment 
begins to oscillate at the same time that the dipole mode 
begins to grow. In the cases where ad = 0.05 or 0.1 when 
the dipole reaches its asymptotic value, the octupole mo¬ 
ment quickly assumes the opposite sign so that the late time 
octupole moment saturates at the same time as the dipole 
but with opposite polarity. This suggests that the lower or¬ 
der multipole moments are excited first, and higher order 
multipoles will be excited and grow as the lower order poles 
reach a saturation value. The octupole moment saturates to 
a value of /x oc t « —N d as suggested by the linear theory. 
For weaker dynamos, the higher order modes are not ex¬ 
cited. These results are consistent with the symmetries of 
the seed field - the disc dynamo is a symmetric function of 
position Z, therefore we expect to excite symmetric (dipole 
and octupole) and not asymmetric (quadrupole) modes. Re¬ 
versing the polarity of the seed field reverses the signs of the 
generated multipoles, but does not otherwise change the re¬ 
sults, as expected. 

We term this the weak dynamo regime because the 
physics is largely determined by what is happening in the 
inner disc. In particular, we performed simulations where we 
varied the cutoff radius 5 ^ i?d yn ^ 8 but kept the dynamo 
coefficient A in ( |23[ ) constant at Rd yn by varying a d . The 
late time moments changed by roughly 10% and growth 
rates and times were similarly affected. Poloidal field plots 
retained the same qualitative character, where a large loop 
formed near the cutoff radius but field lines in the outer 
disc tilted more towards the equatorial plane of the disc. 


3.1.2 Growth Rate 

The disc averaged equations predict that the exponential 
growth rate a oc fiad^K, scaling like yWd and the local 
Keplerian angular velocity. We plot vertically integrated Bz 
at radii R = 10,15 and 20 as a function of angular phase 
flxt for the case ad = 0.1 (Fig. [f^. We scale the magnetic 
field by afi^ 2 and the phase by afi 2 for ease in comparison 
of the other simulations. We see that there is indeed a period 



Figure 3. Dipole moment /tdip as a function of time for dynamo 

numbers Nd = 1 (red, solid), 5 (black, dot) and 10 (green, dash). 

These runs had a constant disc diffusivity a v = 0.1 and the 

dynamo number was varied using the «d y n parameter. By scaling 

1/2 — 1/2 
the time axis to a. d and the dipole moment //dip to a d we 

see the universal scaling of the late time dipole moment and of 

the growth rate. 



Figure 4. Quadrupole moment /i qU ad as a function of time for 
dynamo numbers Nd = 1 (red, solid), 5 (black, dot), 10 (green, 
dash). These runs had a constant disc diffusivity =0.1 and 
the dynamo number was varied using the adyn parameter. We 

have scaled the time axis to a d and the quadrupole moment 

. - 1/2 
/Iquad f° a d 


of exponential growth lasting a time scale of 0(1) where the 
field likewise changes by 0(1). The scaling ( |15| ) suggests 
that a « \/'i in these units but we find, averaging the rates 
for these three radii and ad = 0.01,0.05 and 0.1 runs a = 
2.0 ± 1.0. This suggests that our linear analysis where we 
dropped order one numbers is valid during this short time. 
At late times a gap begins to form for R < Rd yn , as the 
inner disc becomes magnetically dominated owing to the 
advected magnetic flux. 

The magnetic field strength in the disc is also found to 
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increase in time. We define the averaged field in the disc via 


Bi = 


'P>Pflo 


{Bi) 2 dV 


1 P>P flo 


dV 


1/2 


(29) 


where pn 00 r = 1. We define the averaged field in this way 
because Br and B$ are antisymmetric about the disc plane 
and average out to zero vertically. We show a plot of Br/B$ 
and Br as functions of time in Fig. [7]and Fig. [5]. We see that 
the radial magnetic field is an increasing function of time 
and grows nearly linearly. We see that the ratio of radial to 
toroidal magnetic field is, to an order one number, constant 
over the simulation time. This suggests that the simple 
picture whereby differential rotation creates toroidal field, 
and the a-dynamo transforms this field into poloidal field 


is qualitatively correct. Unlike on the surface of the star, 
where the field saturates at late times in the simulation, 
the integrated field in the disc continues to grow. Since the 
growth rate is set by the local Keplerian time, on time scales 
of the simulation the outer part of the disc does not have 
time to saturate with magnetic field, whereas near the star, 
where time scales are expected to be shorter and field is ad¬ 
verted inwards at a faster rate the field has time to saturate. 
We can estimate the growth rate of B by assuming a local 
exponential growth rate at^ln and initially uniform field 
Bo in which case f jj 2 e 2at ^ GMR / RdR ~ l 4 / 3 . If the field 
were decreasing radially then this dependence on t would be 
weaker. The greatest contribution of this integral occurs for 
aflt > 1 so there is a characteristic radius f? max oc be¬ 
yond which there is little contribution. We find B ~ l 11 so a 



Figure 5. Octupole moment /toct as a function of time for dy¬ 
namo numbers Nj = f (red, solid), 5 (black, dot) and 10 (green, 
dash). These runs had a constant disc diffusivity o: r/ = 0.1 and 
the dynamo number was varied using the o: ( j parameter. Note 
how the octupole component has the opposite sign of the dipole 
component at late times. 



Figure 6. Vertically integrated magnetic field Bz at radii R = 
lO(solid), 15(dash) and 20(dash — dot) as a function of angular 
1 /2 

phase Qt scaled by a d f° r the case a^ = 0.1. The exponential 
growth predicted by the linearized regime occurs for phase angles 
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Figure 7. Average radial magnetic field strength Br computed 
from |29| for the cases ad = 0.01 (red), 0.05 (green) and 0.1 
(black). Unlike the field around the star, the disc averaged mag¬ 
netic field does not have time to saturate on time scales of the 
simulation. 



Figure 8. Average radial magnetic field strength Br/B^ com¬ 
puted from [29] for the cases o,j. = 0.01 (red), 0.05 (green) and 
0.1 (black). This ratio varies by an order one number, suggesting 
that the picture where toroidal field is built up via differential 
rotation and a fraction of this is converted to poloidal field via 
dynamo action is qualitatively correct. 


slightly stronger dependence than predicted by this model. 
However, it provides a semi-quantitative understanding of 
how polynomial growth is achieved globally whereas growth 
is exponential locally. 


3.1.3 Outflows 

An interesting question is how the field generated by the dy¬ 
namo affects matter outflows. We measure the total mass 
flux M at the outer boundaries Z = ±5 for various dynamo 
numbers (Fig. [9|. Choosing to measure mass flux along 
boundaries in the range 3 < Z < 8 had little effect on the 
results. We also require that any matter flux have a velocity 
|i> p | > 0.1, which corresponds to an observationally interest¬ 
ing v ~ 200km/s. We note that non trivial outflows are gen¬ 


erated for dynamo number Nd = 10. The outflows begin at 
approximately t = 150 inner disc orbits, corresponding to a 
time t ~ 5yrs for T Tauri stars. This is when the dipole mo¬ 
ment reaches its saturation value and the plasma fl < 1 in 
the inner parts of the disc near R = Rd yn ■ Disc dynamos are 
therefore suitable candidates to explain time varying out¬ 
flows observed in T Tauri systems with observed time scales 
between outflows of 10 — lOOyrs. Such a model has previ¬ 
ously been explored by Stepanovs, Fendt and Sheikhnezami 
(2014) where the dynamo was explicitly turned on and off 
to produce time varying outflows. Though oscillatory, the 
mass flux M « 0.025 or approximately 6.8 x 10 _9 MQ/yr 
for a T Tauri star, which is consistent with observational 
values. The start of outflow generation corresponds to a 
decrease in the ratio of matter to magnetic pressure ft. We 
plot this parameter, averaged over the thickness of the disc, 
as a function of radius for t — 0, 100, 200 and 300 in Fig. 
|10| Winds are primarily launched from the inner part of 
the disc. We see that the region R, < R < 5 becomes mag¬ 
netically dominated (/3 < 1) for t > 200, consistent with 
the formation of a gap for R < Rd yn - The time of outflow 
generation also corresponds to the time at which the stel¬ 
lar multipoles asymptote to their late time values. At this 
time our linear analysis (101 breaks down because we can no 
longer ignore the vz velocity component. Flux can now flow 
in the Z direction, thereby halting the growth of flux on the 
surface of the star. We note that though the stellar dipole 
and octupole moments saturate, other features of the disc 
remain non-stationary - the stellar mass flux and mass flux 
in the wind continue to vary with time. This is consistent 
with the picture that winds are primarily affected by the 
disc, in particular near the inner disc where the mass flux is 
greatest. Outflows can be enhanced via dynamo growth on 
time scales of the simulation. Cowling’s theorem can be cir¬ 
cumvented because the timescales for magnetic field growth 
are shorter than the diffusion timescale. 


3.1.4 Dipole Seed Field 

The time when significant outflows are generated corre¬ 
sponds to the time when the dipole moment reaches its 
saturation value of /Xdi P « 3. We compare this to a run 
where the initial magnetic field is a dipole with magnetic 
moment /Xdi p = 3, corresponding to this late time magnetic 
moment of the dynamo generated field. We keep other pa¬ 
rameters unchanged, ad = a, = a„ = 0.1. This field config¬ 
uration generates outflows roughly an order of magnitude 
smaller than the late time dynamo configuration, but at 
earlier times - since the stellar dipole moment is already 
present, it can begin generating outflows right away with¬ 
out first having to be generated in the disc. Torkelsson and 
Brandenburg (1994a) showed that if a stellar dipole is used 
as the seed field, the field generated in the disc has opposite 
polarity. Dyda et al. (2015) showed that anti-aligned stellar 
dipole and disc fields had suppressed outflows compared to 
a purely disc field because the field line structure in the in¬ 
ner part of the disc was not conducive to mass loading. Our 
results here are consistent with this, where the case with a 
seed disc field has the strongest outflows. If the dynamo is 
turned off, this stellar dipole configuration settles to a con¬ 
figuration with /Tdip ~ 2.5 and the higher order moments 
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Figure 9. Total mass flux A I as a func¬ 
tion of time for dynamo numbers Ay/ = 

0 (black, dotted), 1 (black, solid), 5 (green, dash — dot), 10 (red, dash) 
{ad = 0, 0.01, 0.05, 0.1). Observationally significant outflows 
occur after t = 150. 



Figure 10. Plasma /3 = 87tP/|B| 2 as a function of radius for t 
= 0 (solid), 100 (dash), 200 (dash-dot) and 300 (dotted) in the 
case ad = 0.1. Note that when the inner disc R < 5 becomes 
magnetically dominated, an inner disc wind begins to outflow. 


are zero. If the dynamo is turned on, the dipole moment 
grows and saturates to pdip ~ 4, and the octupole moment 
saturates to p oc t ~ —0.5. We observe the same behaviour 
as for the disc field where the dipole and octupole have the 
same sign until saturation, at which time it quickly reverses 
direction. Though the disc dynamo generates a dipole com¬ 
ponent of the magnetic field, it is unclear whether a disc dy¬ 
namo can generate a closed magnetosphere which has been 
shown to generate strong outflows in the propelor regime. 


3.1.5 Magnetic Diffusivity 

We vary the dynamo number by keeping ctd = 0.1 fixed and 
varying the magnetic diffusivity from 0.03 ^ a v ^ 0.3. The 
viscosity is also kept constant at a„ = 0.3. Below we plot 
the dipole moment as a function of time for these various 
runs. As before, the dipole moment saturates to a certain 



Figure 11. Dipole moment Pdip as a function of time for fixed 
ad = 0.1 and viscosity a v = 0.3 with varying magnetic diffusivity 
a, j = 0.1 (black, solid), 0.07(red,dash), 0.05 (green, dash-dot) 
and 0.03 (blue, dotted). 


value. Interestingly, the sign of the dipole moment depends 
on the magnetic diffusivity. For a v ^ 0.05 then p dip > 0 
and for a v *2 0.07 then pdi P < 0 (Fig. 111. Empirically the 
field line structure is different in these two regimes. In the 
case of larger diffusivity, the field lines in the outer part of 
the disc are mostly vertical, threading through the disc. In 
contrast, for the smaller diffusivity cases the field lines are 
nearly parallel with the disc. This is what we would expect 
as field lines will naturally tilt towards the disc if they are 
frozen in to the matter and the matter accretes inwards. 
Decreasing the dynamo parameter ad = 0.01 had no effect 
on the sign of the generated dipole moment - in fact, the late 
time value pd oc a d holds for ad < 0.1 suggesting that 
the effect is not due to a changing dynamo number. Rather, 
the sign of the late time dipole moment depends on the sign 
of the poloidal field at early times near the star. This is 
controlled by the relative strengths of viscous and diffusive 
effects. The sign change in the generated dipole is likely due 
to an increase in viscous effects in the disc. The Prandtl 
number in the case where a large positive dipole grew was 
V = 10 whereas for other runs this ratio was V < 1. In cases 
where the late time dipole moment asymptotes to a negative 
value, the poloidal field in the upper hemisphere is negative 
at early times. This is reversed for cases where the dipole 
moment asymptotes to a positive value. As described in 
Section [2.2| changing the relative strength of the viscosity 
and diffusivity changes the sign of a term governing the 
flux function evolution. Since the flux function is coupled 
to the toroidal field, we expect this to change the sign of 
the toroidal field. 


3.2 Strong Dynamo 

For ad > 0.1 the evolution of the disc is qualitatively differ¬ 
ent. We set a v = «„ = 0.1 and take ad = 0.2,0.28 and 0.5. 
Below we show poloidal plots of the disc density p, mag¬ 
netic flux contours T and poloidal velocity vectors v p of a 
typical run in the strong dynamo regime in Fig. [T2] for the 
representative case of ad = 0.28. Magnetic loops form for 
all radii R > Rdyn unlike in the weak dynamo regime where 
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loops only form close to -Rdyn- Though loops form along the 
entire length of the disc, the magnetic pressure is strongest 
in the inner portions of the disc. This innermost loop opens 
up due to the differential rotation of the disc and pushes out 
the flux from the outer disc. Because the dynamo is much 
stronger, the disc becomes magnetically dominated and un¬ 
stable. The dynamo time scale is too short compared to the 
viscous and diffusive time scales and disrupts the global 
evolution of the disc on orbital time scales. 

Suppose that to become dynamically significant, the 
seed field must grow by some factor / > 1. Then from (15) 
we get rapid dipole mode growth in one rotation period if 
a.d{3 + ^oid) > (//2tt) 2 . To become magnetically dominated 
/ ~ — | ln/3(0) ~ \ ln(10 3 ) « 3.5 using typical simulation 
seed field values. This yields ad > 0.09, so roughly what 
was found empirically for the transition from the weak to 
the strong dynamo regime. 

As the field grows and dominates the disc dynamics 
on orbital time scales the disc warps and breaks up as a 
gap opens near R ~ Rdyn where the magnetic pressure is 
greatest. Because the disc becomes unstable and breaks up 
we can only study these cases for t ~ 100 inner disc orbits. 
The linear theory, (141 suggests that the growth rate croc aj 
when ad ~ 1. However this is not seen in our simulations. 
In this regime we do not see a period of exponential growth 
suggesting that if it occurs, it is for such a short period of 
time as to be unnoticeable. The gap can open on much faster 
time scales owing to an accretion rate of order a few larger 
than in the weak dynamo regime. The loop like structure of 
the generated field causes the field components to alternate 
signs in the disc plane, and field anihillations lead to chaotic 
behaviour. 


4 CONCLUSION 

We investigated the effects of an a-dynamo operating 
within an accretion disc. Simulation results were affected 
by three dimensionless parameters controlling the dynamo 
strength ad, the disc viscosity a v and the magnetic diffu- 
sivity a v - Our main findings are summarized below. 

1. An initial poloidal seed magnetic field threading the disc 
grows via the disc dynamo - the differential rotation of the 
disc creates a toroidal magnetic field, which is converted to 
poloidal magnetic field via dynamo action. This process re¬ 
peats until at late times the poloidal field strength saturates 
to a value oc yfad. The local growth rate of the field during 

early times is exponential and set by the local Keplerian 

1 / 2 

velocity and the strength of the dynamo with a oc SIk a d ' 

2. The dynamo excites dipole and octupole modes, allow¬ 
ing for the growth of these modes on the star. Quadrupole 
modes are suppressed at late times, owing to the symmetry 
of the initial seed field. Dipole and octupole modes satu¬ 
rate to values of opposite polarity with magnitudes oc fiod. 
This has interesting observational consequences since a disc 
dynamo could be responsible for the observed dipole and 
higher order multipoles observed on T Tauri stars. 

3. A seed field that is too small to generate observationally 
interesting outflows can grow on time scales r ~ 

to values which are high enough to launch observationally 
significant outflows. This suggests that disc dynamos may 


be important in generating the large scale, ordered mag¬ 
netic fields needed for MHD launched outflows. 

4. A stellar dipole seed field generates a field with outflows 
~ 10 weaker than a disc seed field because of interactions 
between the late time disc field and the stellar dipole. This 
may be useful in distinguishing between systems that ad- 
vected their magnetic field from the ISM and those that 
generated them local via a dynamo mechanism. 

5. Increasing the Prandtl number causes the polarity of the 
late time dipole moment to change. The disc dynamo phase 
space is more complex than simply ad as might naively be 
thought - it is jointly a function of three time scales: the 
dynamo, viscous and diffusive times. 

6. For larger values of the dynamo parameter ad > 0.1 
the disc becomes magnetically dominated on time scales 
t ~ f - a gap forms in the inner disc where the dynamo 
is strongest, eventually leading to the break up of the entire 
disc as the growth of the field goes unchecked. 

Physically these findings can be reconciled by under¬ 
standing that three time scales are at play - the dynamo 
time scale, the viscous time scale and the diffusive time 
scale. If the dynamo time scale dominates, the viscous disc 
becomes magnetically dominated and breaks up. This is 
the “strong dynamo” regime and corresponds to the states 
described in Stepinski & Levy (1988) who concluded that 
strong local dynamo modes could disrupt the global modes 
of the disc. If the diffusive time scale dominates then the 
field does not grow sufficiently to launch outflows. When 
time scales are comparable we are in the “weak dynamo” 
regime where outflows can be enhanced but the overall disc 
dynamics are not disrupted. The relative strength of the 
viscous and diffusive effects determines the sign of the late 
time dipole field. 
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APPENDIX A: NUMERICALLY EXTRACTING 
MULTIPOLE MOMENTS 

Any magnetic field can be described by its multipole ex¬ 
pansion. Below, we derive formulae for the axysymmetric 
dipole, quadrupole and octupole components of a magnetic 
multipole expansion in terms of surface integrals in cylindri¬ 
cal coordinates. These are used to numerically extract the 
multipole coefficients of the magnetic field configurations in 
our simulations. 
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t = 70 


t = 100 


z 




Figure 12. Density p (color), poloidal field lines ’3/ (white) and poloidal velocity v p for the case Nd = 28 = 0.28). t = 0 shows 

the initial large scale seed field. By t = 40 magnetic loops have formed along the entire length of the disc. At t = 70 we can see the 
quadrupolar nature of the field as loops are advected inwards. At late times these loops expand into open field lines. The magnetic 
pressure has grown enough to seriously disrupt the disc and open a gap at Rd yn = 5. 


Al Dipole Moment 


general the dipole moment is given by 


Consider a dipole moment, symmetric about the Z-axis with 
flux function given by 


'3/ — /^dip _ 


R 2 


' ( R 2 + Z 2 ) 3/2 ’ 
and poloidal magnetic fields 

o _o RZ 

-LsR O/idi] 


(7? 2 + Z 2 ) s/2 ’ 


(Al) 

(A2) 


^JdV [r x j] = ^JdW [r x (V x B)], (A4) 


where we work in CGS units and have used Ampere’s law. 
We will integrate over a cylinder 0 < R < a, —h < Z < 
h, 0 < (j> < 2n enclosing the star. We also assume the 
boundary condition i&(0, Z) = 0 Using vector identities we 
may write 


Bz — /r dip 


( 2 Z 2 - R 2 ) 
(R 2 + ^ 2 ) 5/2 ’ 


(A3) 


Given a general magnetic field, we want to extract its dipole 
moment from its multipole expansion at any given time. In 


m = J rfv [v (r ■ B) - (r ■ V) B - (B ■ V) r]. (A5) 

We are interested in the Z component of the magnetic mo¬ 
ment, Udip- We calculate the contribution to fi, di p from each 
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Figure Al. Flux lines of axially symmetric dipole (left), quadrupole (center) and octupole (right). 


of the above three terms below. 

T, = [ dV V z (r ■ B) = [ dS z (r ■ B) 


= 2 ^ r 

Jo 


(A 6 ) 


RdR (r ■ B) 


where we have used the divergence theorem to convert the 
volume integral to a surface integral over the sides of the 
cylinder and only the top and bottom contribute to the Z- 
component. Using the definition of magnetic field in terms 
of the flux function we may write 

Ti = 2 tt^ dRR 2 B R ^+2TTh^H(a,h) + 'l>(a,-h)^. (A7) 
The Z-component contribution from the second term is 
T 2 = - f dV [(r-V)B] 




(A 8 ) 


Integrating the first term by parts and using the definition 
of the magnetic flux function we may write 

T 2 = — 27t J dZ^a 2 Bz(a, Z) — 2'i>(a, Z) — RZB R (a, Z)^. 

(A9) 

The third term can be dealt with by rewriting it in terms of 
the flux function and explicitly carrying out the integration 
in R to yield 


T 3 = -2n [ dZ^(a,Z). 

J-h 

Combining these terms we arrive at the result 


(A10) 


dip — ‘ 


-j^ dRR 2 B H (R,Z)+h^{a,h) + V(a, -ft)] 
- J dZ | a 2 Bz{a, Z) — ^(a, Z) — aZB R (a, Z )j 


which is easy to implement in our cylindrical code. 
© 2015 RAS, MNRAS 000,[lj{l5| 
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A2 Quadrupole Moment Calculation 

Consider a quadrupole moment, symmetric about the Z- 
axis with flux function given by 

R 2 Z 


$ = 


3/Tquad 


4 (R 2 + Z 2 f 2 ’ 

and poloidal magnetic fields 

3/iquad R{R 2 - 4 Z 2 ) 


B r = - 


B z = 


4 ( R 2 + Z 2 ) 7/2 ’ 

3/iquad ( Z 2 — 3 R 2 ) 


(A12) 


(A13) 


(A14) 


4 (R 2 + Z 2 ) r/2 ' 

Its quadrupole moments are Du = D 22 = — D 33/2 where 


/^quad — D 33 — 


/dVt 


r • (V x J)] Z 2 


(A15) 


= 2 / dV Z (r x J) z 


where as with the dipole, we are interested in the Z- 
component of /U qu ad- Using Maxwell’s equations we can 
write 

fiquad = ^ j dVZ [r x (V x B) ] (A16) 

As with the dipole case, this expression can be expanded 
using a vector identity as 

Mquad = ^ J dVz[v(r-B)-(r-V)B-(B-V)r] z . 

(A17) 

As in the dipole case we calculate each of these terms 
separately. We are interested in the Z component of the 
quadrupole moment. 


Ti = J dV Z V z (r • B) 

= - [ dV (r • B) 2 + Z (r • B) 


(A18) 


where we have integrated by parts. Writing out the mag¬ 
netic field in term of the flux function, we expand the first 
term and carry out some of the intergrals. The second term 
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we evaluate on the top and bottom of the cylinder which 
yields 


Ti = 2 tt 



dR R *{R,z) 


h 

h 


-2tt f dZZV(a,Z) 

J-h 


+ 2n f dR R Z ( RBr + ZBz ) 

Jo 


h 

— h 


(A19) 


Expanding the second term and carrying out the angular 
integral we find 


T 2 = - J dVZ( r ■ V)B 


(A20) 


Performing integration by parts for the R integral on the 
first term and then writing it out in terms of the flux func¬ 
tion allows the integral in R to be evaluated, yielding 


T 2 = 2tv [ dZ 2Z 4>(a,Z) - 
J-h 


a 2 ZB z (a, Z) 


+ 2-7T 


f 


dZaZ 2 Br(cl, Z). 


(A21) 


For the final term, we take the divergence, expand it in 
terms of the flux function, carry out the integral in R and 
find 

T 3 = - J dVZ {B ■ V)Z = -2 tt J dZ Z ^(o, Z). (A 22 ) 

Combining these terms we arrive at the result 

a h 

/^quad — C J dR\m(R,Z) + RZ(RB r + ZB z )]\ 

I J dz[aZ 2 B R (a,Z)-a 2 ZB z (a,Z)]. 

(A23) 


+ c 


A3 Octupole Moment Calculation 

Consider an octupole moment, symmetric about the Z-axis 
with flux function given by 


Meet R 2 (4Z 2 - R 2 ) 

2 (f?2 + Z 2 ) 7/2 ’ 

and poloidal magnetic fields 

= 5/r oct RZ{AZ 2 - 3 R 2 ) 
2 (/,■- . 


(A24) 


(A25) 


B z = 


Moot ( 8 Z 4 - 24 Z 2 R 2 + 3R 4 ) 


2 (R 2 + Z 2 f ' 2 

The octupole moment is given by 


where 


Moct = 2 —M 30 , 


M 30 = i J dV [r x J] ■ V (r 3 y 3 o) , 


(A26) 

(A27) 

(A28) 


and 

r 3 F 3 o = \ ~!~r 3 P 3 ( cose) = ]\ -Z (2Z 2 - 3R 2 ) , (A29) 

V 47T 4 v 7T 

is a spherical harmonic that we have written in terms of 
the Legendre polynomial P 3 . Expanding the cross product 
and using Amperes law to write out the current in terms 
of the magnetic field and evaluating the azimuthal integral 
we find 


l^oct 


3c 

16 


'/■ 


4 / dRd ztfz* - H) 


7 


- / dRdZR 


(dB R 

0B Z \ 

V dZ 

dR J 


(A30) 


Evaluate each of these terms separately. The first term is 

Ti m f dRdZR 2 f 


= / dRR 2 h 2 B R 


dB R 

dB z \ 

dZ 

dR J 

-7 

dR dZRV + 2 [ dR,RZ4> 

h 

J 

Z) + 2 

J dZZ 2 4>(a,Z), 


(A31) 


where we have used integration by parts and the definition 
of magnetic field in terms of the magnetic flux function. 
Similarly, the second term is 

Es / dJMZR4 (w + w) 

h 

= f dRR 4 B R - f dZa 4 Bz(a,Z) (A32) 

J -h ' 

-8 J dRdZR'& + 4 J dZa 2 ^(a,Z). 

Combining these terms we find 
3c 


/^oct — 


16 


7 


4 / dR R h Br\ +8 dR RZ'i 


!■ 


4 J dZ Z 2 a 2 B z {a, Z) + 8 J dZ Z 2 4>(a, Z) 
- J dR R 4 B r + J dZa 4 B z {a,Z) 


-4 j dZ a 2 'L(a, Z) 


(A33) 


We note the cancellation of the volume integrals in Ti and 
T 2 , allowing the final expression to again only depend on 
surface integrals. 
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